In this paper, an iterative method is constructed to solve the following constrained linear matrix equations system: [A1(X),A2(X),... ,Ar(X)]=[E1,E2, ... ,Er ], X ∈ I={X |X= U(X)}, where Ai is a linear operator from Cmxn onto Cpixqi, Ei ∈ Cpixqi, i=1 , 2,..., r , and U is a linear self-conjugate involution operator. When the above constrained matrix equations system is consistent, for any initial matrix X0 ∈ I, a solution can be obtained by the proposed iterative method in finite iteration steps in the absence of roundoff errors, and the least Frobenius norm solution can be derived when a special kind of initial matrix is chosen. Furthermore, the optimal approximation solution to a given matrix can be derived. Several numerical examples are given to show the efficiency of the presented iterative method. Mathematical subject classification: 15A24, 65D99,65F30.